Repository navigation
Fix fmaf() double-rounding error (#160) - #357
Open
ViralBShah wants to merge 1 commit into
Open
ViralBShah wants to merge 1 commit into
ViralBShah wants to merge 1 commit into
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #357 +/- ##
==========================================
+ Coverage 72.12% 72.14% +0.01%
==========================================
Files 233 233
Lines 6135 6139 +4
Branches 1607 1607
==========================================
+ Hits 4425 4429 +4
Misses 1417 1417
Partials 293 293 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
fmaf(0.9474001f, 4.639901e-7f, -0.24325085f) returned -0x1.f22d44p-3f (-0.2432504f) instead of the correctly single-rounded -0x1.f22d46p-3f (-0.24325041f), a 1-ulp double-rounding error. Root cause: fmaf computes (double)x*y + z and, for halfway-to-float results, applies a correction unless the double result is exact. The exactness check used "result - xy == z". That test is only valid when |xy| >= |z|; when |z| > |xy| the difference result - xy is itself rounded and can equal z even though result already dropped low-order bits of the true sum. The false "exact" verdict skipped the correction and let the inexact double result be rounded a second time to float -> double rounding. Replace it with Knuth's 2Sum, which yields the exact rounding error of xy + z independently of operand magnitude; result is exact iff err == 0. volatile barriers keep the decomposition exact under FP contraction. Verified: the reproducer now returns -0x1.f22d46p-3f, and the rebuilt library matches a hardware (single-rounded) FMA on 150M+ random/edge triples (normal-range and full random bit patterns) with zero mismatches. Full make test passes (test-double, test-float, regression). Adds test/regression/test-160.c pinning the reproducer and several hardware-FMA-verified triples (teeth-checked: fails on the old code). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
r3bb1t
pushed a commit
to r3bb1t/rusty_box
that referenced
this pull request
Oct 26, 2026
f64_fused/f32_fused trusted Rust's mul_add on every std target. Where it lowers to a C library's fma/fmaf, some round normal results wrongly (MinGW-w64 on MSVCRT, rust-lang/rust#140515; FreeBSD msun s_fmaf.c, used by bionic on non-Arm Android ABIs, JuliaMath/openlibm#357), and the result bound does not filter those errors. HOST_FUSES_ONCE, next to HOST_ROUNDS_ONCE, is the allow-list; every other host answers None and SoftFloat computes (R5: one choke point for the hazard). Allow-list and evidence: - target_arch = "aarch64": llvm.fma is FMADD/FMSUB of the base Armv8-A FP, which rounds once (Arm ARM FPMulAdd). - x86/x86_64 with target_feature = "fma": llvm.fma is VFMADD, single rounding (Intel SDM). - x86_64 linux gnu: on this host the perfbench (task6, task6fix) and the lib test binary import no fma/fmaf (nm -D); both are local symbols from compiler-builtins' libm port (arch/x86/fma.rs: runtime dispatch to vfmadd213sd/ss, FMA4 vfmaddsd/ss, else generic fma_round = musl fma.c, correctly rounded for all inputs, and fma_wide_round = musl fmaf.c, correctly rounded wherever the result is normal, which the bound requires). glibc's own fma/fmaf (dbl-64/s_fma.c, flt-32/s_fmaf.c) are correctly rounded for every input. New test the_host_fuses_after_the_forms_negations: with FPSR.IXC set, FMADD/FMSUB/FNMADD/FNMSUB through f64_mul_add/f32_mul_add give SoftFloat's bits and status for the hand-negated operands, and the host answers for those operands exactly where HOST_FUSES_ONCE holds. Mutation: moving the host call above the negations (f64; and above Rn's negation, f32) fails it on Fmsub; no other test caught either mutant. The random fused test's reach thresholds now apply where HOST_FUSES_ONCE holds, else accepted == 0. icount (run_until, callgrind): fp plain 12,400,149,017 and hooked 13,160,148,837, identical to 5e71b4b; hashes identical to d21b7e9. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Jr1yEuXJkJ5uLHvK7KJWJ8
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fixes #160.
fmaf(0.9474001f, 4.639901e-7f, -0.24325085f)returned-0x1.f22d44p-3f; the correctly single-rounded result is-0x1.f22d46p-3f(matches a hardware FMA).Root cause
src/s_fmaf.ccomputes(double)x*y + zthen corrects halfway cases — but usedresult - xy == zas an exactness shortcut. That test is bogus when|z| > |xy|(here|xy|≈2^-22,|z|≈2^-3):result - xyis itself rounded and returnszeven thoughresultalready dropped low-order bits. The false "exact" verdict skips the correction, and the inexact double is rounded a second time to float — a 1-ulp double-rounding error. (FreeBSD/msun upstream has the same bug.)Fix
Replace the shortcut with Knuth's 2Sum, which yields the exact rounding error of
xy + zregardless of operand magnitudes (resultis exact ifferr == 0).volatilebarriers keep each step rounding in double so it stays exact even under FP contraction.Verification
-0x1.f22d46p-3f.-mfma -ffp-contract=fast.make testpasses.test/regression/test-160.c(float; runs on all arches). Fails before / passes after.